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Abstract 

We give an algorithm to compute the series expansion for the inverse of a given 
function. The algorithm is extremely easy to implement and gives the first terms of 
the series. We show several examples of its application in calculating the inverses of 
some special functions. 

MSC-class: 30B10 (Primary) 30D10, 33E20 (Secondary) 



1 Introduction 

"One must always invert." 
Carl G. J. Jacobi 



The existence of series expansions for inverses of analytic functions is a well-known result 
of complex analysis ^7]. The standard inverse function theorem, a proof of which can be 
found, for example, in states that 

Theorem 1.1 Let h{x) he analytic for \x — xq\ < R where /i'(xo) 7^ 0. Then z = h{x) has 
an analytic inverse x = H{z) in some e -neighborhood of zq = h{xo). 

*dominicd@newpaltz.edu 
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In the case when xq = zq = 0, \h{x)\ < M for \x\ < R, and h'{0) = a, R. M. Redheffer 
has shown that it is enough to take e — - ^"'^'^ 



4 M 



However, the procedure to obtain the actual series is usually very difficult to implement 
in practice. Under the conditions of Theorem ll.il the two standard methods to compute the 
coefficients 6„ of 

h~\z) = Hiz) = j2uz-^or 



n>0 



are reversion of series ^H], [2^]; and Lagrange's theorem. The ffist one requires to 
expand h{x) around xq 

h{x) = a„(x — xo)"' 



and then solve for bn in the equation 



n>0 



n>0 



n>0 



Z - Zo) 



Xo 



by equating powers of z and taking into account that oq = zq and bo = xq. This method 
is especially useful if all that is known about h{x) are the ffist few a„. When xo = zq = 
and tti = a, it was shown by E. T. Whittaker |34j that 



a 



, a2 

02 = T7. 



3!a5 



l„2n-l 



n\a 



na2 a 
2na3 {n + l)a2 
3raa4 (2n + l)a3 + 2)a2 



3a2 a 
603 4a2 



2a 






3a 



Ana^ (3n + l)a4 2(n+l)a3 (n + 3)a2 



where |-| = det (■) . In Example 10, we show how to get the 6„ in term of the a„ using our 
method. 

A computer system like Maple can reverse the power series of h{x), provided h{x) is not 
too complicated, by using the command 

> Order := + 1; 

> solve ( series ( h{x) , x = xo , N + 1) = z , x); 

where N is the number of terms wanted. Fast algorithms of order (nlogn)^/^ for reversion 
of series have been analyzed by Brent and Kung jH], [SI- The multivariate case has been 
studied by several authors jlj, [Sj, [E], and Wright jSHI has studied the connection 
between reversion of power series and "rooted trees" . 
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The second and more direct method is Lagrange's inversion formula PP, 



1 d 



n\ dx 



n~l 



X — Xq 

h{x) - zq 



x=xo 



Unfortunately, more direct doesn't necessarily mean easier and, except for some simple cases, 
Lagrange's formula (jl.ip is extremely complicated for practical applications. The g-analog 
(a mathematical expression parametrized by q which generalizes an expression and reduces 
to it in the limit q 1+) of p.ll) has been studied by various authors [2], [H], [IH], 120] 
and a unified approach to both the regular and g-analog formulas have been obtained by 
Krattenthaler There has also been a great deal of attention to the asymptotic expansion 
of inverses IZT], [2H1, EH, E2]- 

In this note, we present a simple, easy to implement method for computing the series 
expansion for the inverse of any function satisfying the conditions of Theorem II. H although 
the method is especially powerful when h{x) has the form 

X 

h{x) = / g{x)dx 



and g{x) is some function simpler than h{x). Since this is the case for many special functions, 
we will present several such examples. This note is organized as follows: 

In section 2 we define a sequence of functions S}"[/] (x), obtained from a given one /(x), 
that we call "nested derivatives" , for reasons which will be clear from the definition. We give 
a computer code for generating the nested derivatives and examples of how ID'^[/] {x) look 
for some elementary functions. Section 3 shows how to compute the nested derivatives by 
using generating functions. We present some examples and compare the results with those 
obtained in Section 1. 

Section 4 contains our main result of the use of nested derivatives to compute power 
series of inverses. We test our result with some known results and we apply the method for 
obtaining expansions for the inverse of the error function, the incomplete Gamma function, 
the sine integral, and other special functions. 

2 Definitions 

Definition 2.1 We define (x), the n*'* nested derivative of the function f{x), by the 

following recursion: 

^V] W = 1 

2)"[/](x) = ^[/(a;)xS)"-i[/](x)], n>l. (2.1) 
Proposition 2.2 The nested derivative 2D'^[/] (x) satisfies the following basic properties. 
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(1) For n > 1, S!)"[k] = 0, k constant. 

(2) Forn>0, S)"[k/] (x) = (x), k constant. 

(3) For n > 1, {x) has the following integral representation: 



•where Ck is a small loop around x in the complex plane. 

Proof. Properties (1) and (2) follow immediately from the definition of (x). 

To prove (3) we use induction on n. For n = 1 the result follows from Cauchy's formula 

Ci 

Assuming that the result is true for n and using ()2.H) . it follows that 

Cn + l 

1 / / /(^^n+l) fiZn) 



(27rz)"+i / / - xy {zn - z^+i)'^ 

X TT ^^^^^^^—-rdZn+l ...dZi 

[Zk - Zk+lY 



k= 



n- — ^ — r^-C^^n+l ■ • • dZi- 



(27rz)"+i J J - x)2 J-^-L [zk - 2:^+1 



Algorithm 2.3 TheD algorithm. The following Maple procedure implements the recurrence 
relation \2.1\) . We define d{k) = T)^[f] (x), where N is the number of terms desired. 

> d{0) := 1; 

> for k from to N do: 

> d{k + 1): = simplify ( diff ( f{x) * d{k), x )): (2.2) 

> print ( k + 1, d{k + 1) ): 

> od: 
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Example 2.4 The function f{x) — x. 



D^[/] (x) = l 
^V] (^) = 1 

^V] (^) = 1- 

Example 2.5 T/ie power function f{x) — x^ , r ^ 1. 

3)1 [/] (a;) = ^2:"-i 

= r(2r - l)a;'("-^) 
D3[/] = r(2r - l)(3r - 2)x^("-') 



= (r - r^x''^'-'-\ 

r(i + ^) 

Notice that when r = -rrr, k = 1,2,..., the sequence of nested derivatives has only 



k+l 

k + 1 non-zero terms 



X fc+i , <n < k 
0, n>k+l 



Example 2.6 T/ie exponenetial function f{x) — e^^ . 

{x) = re''^ 
2)2 [/] (x) = 2r'e2'-^ 

3 Generating functions 

Generating functions provide a valuable method for computing sequences of functions defined 
by an iterative process; we will use them to calculate S)"'[/] (x). In the sequel, we shall 
implicitly assume that the generating function series converges in some small disc around 
z = Q. 
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Theorem 3.1 Given h{x) = J j^dx, its inverse H{x) = h ^{x) and the exponential gen- 
erating function G{x,z) = (x)^, it follows that 



n>0 



Gix,z) = ^ifoH)[z + hix)]. (3.1) 



Proof. Taking ()2.ip into account gives 



- [fix) X Gix, z)] = J2d^ [fix) X ix)] - 

n>0 

n>0 ■ n>l ^ '' 



n>0 

Hence, the generating function (^(x, z) satisfies the PDE 

dif X G) _ dG 
dx dz 

with general solution 

Gix,z) = ^g[z + hix)] (3.2) 

where giz) is an arbitrary analytic function. Invoking the boundary condition Gix^Q) = 
D0[/] (a;) = 1, (j^ gives 

and therefore 

fix) = igoh) ix). 

If X = Hiw), then 

(/ oH)iw) = igoho H) iw) = giw) 
and the theorem follows. ■ 

Example 3.2 The function fix) = x. 

Here /i(x) = J ^dx = ln(x), if(x) = e^, and from 4,'/. i)) it follows that 

Gix, z) = — exp [z + Inix)] = e^. 



We could obtain the same result from Example \2.4\ by summing the series 

z"' 



n>0 
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Example 3.3 The power function f{x) = , r ^ 1. 

r .1 — r 1 

Now h{x) = J x~^dx = ^37, H{x) = [(1 — r)x]'^-'' , and we get 



G{x, z) = X 



'l-r)\ z + 



X 



'1 — r)z + X 



1-r 



X 



1-r 



1 — r 

= [1 + (1 - r)x'-^z\^^ . 

Expanding in series around 2; = 0, we recover the result from Example 

U 1Z7 = ^) i.e. r = A; = 0, 1, . . then G{x, z) is a polynomial of degree k in z 

and hence 

D"[/](x) = 0, n>k + l 
as we have already observed in Example\ 



Given the particular form of the function h{x) in Theorem 2, we can get alternative 
expressions for p.ip which sometimes are easier to employ. 

Corollary 3.4 Let h{x) = f jj^dx, its inverse H{x) = h~^{x) and the exponential gener- 
ating function G{x, z) = (x)^. Then, 

n>0 



(i) 



and 



(ii) 



G{x,z) = -—H' [z + h{x)] 
fix) 



G(x,z) = —H[z + h(x) 
dx 



(3.3) 



Proof. 

(i) By definition {h o H){x) = x, so 
Since h{x) = J j^dt, 



h! [H{x)] H\x) = 1. 



-H'ix) = 1 



or 



and therefore 



flHix)]- 
UoH){x) = H'{x) 

G{x, z) = -^{f oH)[z + h{x)] 



fix 
1 



■H' [z + h{x)] . 
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(ii) Using the chain rule 

-^H [z + h{x)] =H'[z + h{x)] h'{x) 
ax 

= I, + ,(,)]_!. 

and the conclusion follows from part (i). 



4 Applications 

We now state our main result. 

X 

Theorem 4.1 Let h{x) = J jjf^dt, with f (a) 7^ 0,±cxd, audits inverse H{x) = h^^{x). Then, 

a 

~n 

H{z) = a + f{a)J2^-'[f]{a)- (4.1) 
where \z\ < e, for some e > 0. 

Proof. Let's first observe that since h{a) = 0, it follows that H{0) = a and from ()3.3|) 

G{a, z) = j^^H' [z + h{a)] = j^^H'{z) 



where 



Hence, 



G{a,z) = J2^-[f] ia)-j. 



n>0 



H{z) = H{0) + I /(a) (a)^rft 

= a + /(a)5^2)"[/] (a)^^^ 
+ /(a) {c 



n>0 

la)— r 
n! 

n>l 
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Example 4.2 The natural logarithm function. Let f{x) = e ^, with a = 0. We have /(O) 
1, 



h{x) = I e'dt = - 1, H{x) = ln(x + 1) 

and from Example[ 



S)"[/] (0) = {-irnl 
Hence, from we get the familiar formula 



^-^ n 

n>l 

Example 4.3 The tangent function. Let f{x) = x"^ + 1, with a = O.Now /(O) = 1, 

X 

f 1 

h{x) = / -dt = arctan(x), H{x) = tan(2;) 



and d^.il ) implies that 

tan(z) =^2)"-^[x2 + l] (0) — . 
n>i ^■ 

Therefore, 

2)2'=+i[x2 + l](0) = 0, A;>0 
D2^[x2 + l](0) = -^4^(4^+i-l)|52(fc+i)|, A;>1 (4.2) 

where Bj. are the Bernoulli numbers 

Remark 4.4 From Example \2. 5L we recall that 

{x) = + 

and consequently 

(0) = 0, n > 1. (4.3) 

Comparing (|^.ij| ) and (|^.5| ) we can see i/ie highly nonlinear behavior of the nested deriva- 
tives, since even the addition of 1 to f{x) creates a completely different sequence of values, 
far more complex than the original. 

We now start testing our result on some classical functions. 



9 



Example 4.5 Elliptic functions. Let f{x) = a/1 — sin^(x), < p < 1, a = 0. We 
have, /(O) = 1 and 



de 



1 -p2 sin^(^) 



F(p; 0), -f^(p; = am(p; x) 



where F{p; cj)) is the incomplete elliptic integral of the first kind, and am{p; x) is the elliptic 
amplitude 



am(p; x) = arcsin [sn{p; x)] = arccos [cn{p; x)] = arcsin 



-y/l — dn'^{p; x) 
p 



with sn{p; x) , cn{p;x), anddn{p;x) denoting the Jacobian elliptic functions. 
Computing (0) with \2.^) gives 



S^^^M/] (0) = 0, k>Q 

^''^[f] (0) = (-i)ygfcb), k > 1 

where Qk{p) is a polynomial of degree 2{k — 1) of the form 

Q,(p)=p2{^-i) + ... + 22('=-i). 

The first few Qkip) are 

Qi{p) = 1 

Q2{P) =P^+4: 

Qsb) =P^ + 4V + 16 

Q^(p) = / + 408/ + 912p2 + 64 

Qr^(p) = / + 3688p^ + 307682/ + 15808/ + 256 



an 



d implies that 



am{p; x) = z - /— +/(/ + 4)— - /(/ + 44/ + 16)— + 
3! 5! 7! 

in agreement with the known expansions for am{p; x) |,jJJ /. 



(4.4) 



Example 4.6 The Lambert- W function. Let f{x) = e ^{x + 1) ^, a = 0, f{a) = 1. 
Here 

h{x) = xe\ H{x) = LW{x) 



10 



where by LW{x) we denote the Lambert- W function fTJ^ . W^ . In this case, gives 



2)1 [/] (0) = -2 
(0) = 9 
(0) = -64 

(0) = [-(n+l)r. 

From we conclude that 

n 

LW{z) = ^(-l)"-ir2"-i^. 

n>l ^' 

Example 4.7 We now derive a well known result 11] about reversion of series. If we take 

h{x) = aix + a2X^ + a^x^ + a^x^ + a^x^ + agx^ + ayx'' + • • • 
where ai 7^ 0, then 

f(x) = = - 

h'{x) ai + 2a2X + Sa^x'^ + Aa^x^ + 5050;'^ + Ga^x^ + Tajx^ + ■ ■ ■ 

a = 0, /(O) = ^ and from \2.i^) we get 
DM/](0) = -2 



1) ^[/](0) = 6^i^Hl!_f^ 

(oi) 

2) 3^/] (0) = 24^^^!^ ~ ^^^^'^ ~ ^'^^2'* 



(ai) 

6(01)^0203 + 3 (0103)^ + 14(02)^ - (01)^05 - 2101(02)^03 



S'[/](0) = 120 

Oi 



Hence, 



1 ^2 2 , 2(02)^-0103 3 6010203 - (Oi)^ O4 - 5(02)^ 4 

II[z) = — z ^z H 5 z H ^ z 

ai (ai) (oi) (oi) 

6(01)^0203 + 3 (0103)^ + 14(02)^ - (01)^05 - 2101(02)^03 5 

/ x9 Z + ■ ■ • . 

Oi 
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Remark 4.8 An explicit formula for the n}^ term is given in Morse and Feshbach \24\ Part 
1 pp. 411-413], 

^ _ 1 , ^y+t+u+- n{n + l)---{n-l + s + t + u^ ) /a2 V /asV 

niaiY s\t\u\--- \CLij \ai) 

s,t,u,... ^ ' ^ ' 

s + 2t + 3m H = n - 1. 

Example 4.9 The Error Function, erf(a;). We now have 

X 

2 /" _ 2 
h{x) = erf(a;) = —= / e~* dt 

VTT J 



/(a^) = , a = 0, /(a) = ^ 

f 0, n = 2k + l, k>0 



and / li^.i^) g'if es 



where 



Ao = 1, Ai = 2, A2 = 28, A3 = 1016, A4 = 69904 
^5 = 7796768, Aq = 1282366912, A = 291885678464, 



From \4-l\ ) 



n>0 ^ ^ 



(2n+ 1)!' 

which agrees with other authors calculations previously published l^, fT^ . fTS^ . UI^ . 
We will now extend ()4.1|) to a more general result. 

X 

Corollary 4.10 Let h{x) = J jjj^dt, zq = h{b), with f{b) 7^ 0, ±00 and its inverse H{x) = 

a 

h^^{x). Then, 

H{z) = b + f{b) J2 (^)^^^T^ (4-5) 

n>l 

where \z — zo\ < e, for some e > 0. 
Proof. We consider the function 

u{x) = h{x) — zq 
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which satisfies u{b) = 0, and its inverse U{x) = u ^(x). Since f{b) 7^ 0, ±00, we can apply 
fl4.1|) to u{x) and conclude that 



n>l 



All that is left is to see the relation between U{z) and H{z). 
Suppose that u{x) = y. Then 

y = u{x) = h{x) — Zq 
h{x) = y + zo 
X = H{y + 2:0) 

and therefore 

U{y)=H{y + zo) 

or 

H{z) = U{z - zo) 

and (j4.5p follows. ■ 

We will now use our results to get some power series expansions that have not been 
studied before. 

Example 4.11 The incomplete Gamma function, 7(z/;a;). We have 

X 



/(z/;x)=eV-^ a = 0. 

Since 

/(^;0) = 



0, < z/ < 1 
00, 1/ > 1 



we can't apply \4-^^ - Choosing 6=1, Zoi^iy) = 7(z/; 1), f{y]h) = e, we conclude from 
that 

i7(.;.) = l + e$:D-M/](l)J'"'°^"^]" 



n>l 



We use to compute the first few !)"[/] (1) and obtain 

(1) = e-QM 
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where Qn{i^) is a polynomial of degree n 



and we can write 



2-v 

36 - 53z/ + 29z/2 - 

245 - 474z/ + 375z/^ - U6u^ + 24z/^ 

2076 - 4967z^ + 5104//^ - 2847z^^ + 874^^ - 120u^ 



H{u; z) = l + ^e"Q„_i(//)i 



n>l 



Example 4.12 The sine integral function, Si(x). In this case 

X 

sin(f) 



h{x) = Si{x) 



-dt, f{x) 



X 



a = 0. 



For this example f{a) is well defined, but to simplify the calculations we choose b 
f , zo = Si(f ) ~ 1.370762. Then, 

m = I, (f ) = Qn(vr) 

where Qn{x) is once again a polynomial 



Qi{x) 


= 1 




Q2{X) 


= 1 + 


1 2 

r 


Q3{X) 


= 1 + 


7 2 
-X 

4 


Qiix) 


= 1 + 


8x^ + —x^ 
16 


Q5{X) 


= 1 + 


61 159 4 

— X H X 

2 16 


Qeix) 


= 1 + 


423 . 1671 

x^ H 

4 16 



an 



d from \4-5^ we obtain 



H{z) 



n>l 



4 I 225 g 

X H X . 

64 



{z - ZoY 
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Example 4.13 The logarithm integral function, \i{x). From the definition 

X 

h{x)^]i{x) = j 1^^^' f{x)^\n{x), a = 0. 



In this case f{a) — — oo, so we must choose b. A natural candidate is b — e, which gives 

f(b) = 1, zq^ li(e) ~ 1.895117816 

with 

Al = l, ^2 = 0, ^3 = -l, ^4 = 2, ^5 = 1 

= -26, Ay = 99, = 90, = -3627 

and we have 

Hiz) = e + Y,A^-^^. 

n>l 

Example 4.14 The incomplete Beta function, B^u, ij,;x). By definition 

X 



and hence 

f{u,n;x)^x^-"'{l-xy-'', a = 0. 
To avoid the possible singularities at x — and x — 1 we consider b — and therefore 

f(u,fi;b)^^2^+^, Zo{u,fi)^B{u,fi;^). 



The D algorithm now gives 



with Qn{^, y) a multivariate polynomial of degree n 

(5l(z/,/x) = jJi-P 

Q2{iy, n) = -2 + u - 4:un + + 2f/ + 2u'^ 

Q^(u, ij,) = {n- - I2vn + 7H-12 + 6z/^ + 7u) 

Qi{v, /x) = 16 - 46z^//^ - 46j/V - - 22/x + 154i/// - 96i//x^ - 96i/V+ 

144z/V^ - 22i/ - 63z/^ + 24i/^ + 46zy^ + 24//^ + 46//^ 
Q5(z/, n) = {n- z/)(120//^ + 326^1^ - 480I///3 + 720i/V^ - 323/^^ - 326z//x^ 

- 362/^ - 480Z/V + 1154i/// - 326z/V - 3231/^ + 240 + 120i/^ 

+ 326i/^ -362i/). 
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and we have 

H{z) = \ + \r^^Y.'^<^^''''^Qr,{u,^,)^-^^^. 

n>l 

Conclusion 4.15 We have presented a simple method for computing the series expansion 

for the inverses of functions and given a Maple procedure to generate the coefficients in these 
expansions. We showed several examples of the method applied to elementary and special 
functions, and stated the first few terms of the series in each case. 

Acknowledgement 4.16 / wish to express my gratitude to the referees for their extremely 
valuable suggestions and comments on previous versions of this work. 
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